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A unique signature of tiie modified Newtonian dynamics (MOND) paradigm is its peculiar be- 
liavior in tfie vicinity of the points where the total Newtonian acceleration exactly cancels. In the 
Solar System, these are the saddle points of the gravitational potential near the planets. Typically, 
such points are embedded into low-acceleration bubbles where modified gravity theories a la MOND 
predict significant deviations from Newton's laws. As has been pointed out recently, the Earth-Sun 
bubble may be visited by the LISA Pathfinder spacecraft in the near future, providing a unique 
occasion to put these theories to a direct test. In this work, wo present a high-precision model of the 
Solar System's gravitational potential to determine accurate positions and motions of these saddle 
points and study the predicted dynamical anomalies within the framework of quasi-linear MOND. 
Considering the expected sensitivity of the LISA Pathfinder probe, we argue that interpolation func- 
tions which exhibit a "faster" transition between the two dynamical regimes have a good chance 
of surviving a null result. An example of such a function is the QMOND analog of the so-called 
simple interpolating function which agrees well with much of the extragalactic phenomenology. We 
have also discovered that several of Saturn's outermost satellites periodically intersect the Saturn- 
Sun bubble, providing the first example of Solar System objects that regularly undergo the MOND 
regime. 
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I. INTRODUCTION 

Since the beginning of the past century, the Solar Sys- 
tem (SS) has been a major lab for testing the validity 
of general relativity (GR). Ranging from the first evi- 
dence of Mercury's perihelion precession [1] up to the 
recent tests of the Lense-Thirring effect [2, 3], the SS 
has played a key role in our understanding of the laws of 
gravity. All SS tests of GR performed so far, however, 
have tested Einstein's theory in the regime of moderately 
strong gravity, i.e. where the total Newtonian accelera- 
tion exceeds the typical values found at galactic and cos- 
mological scales by several orders of magnitude. At larger 
scales, however, it is well-known that GR fails to account 
for the dynamics of galaxies and galaxy clusters unless a 
massive component of dark matter (DM) [4, 5] is taken 
into account. The current lack of knowledge about the 
composition and behavior of this dark component of the 
Universe, which has recently increased with the study of 
the internal dynamics of dwarf galaxies [6] , together with 
the failing attempts to detect DM particles [7] , motivates 
research on different alternative gravity theories such as 
f{R) gravity [8], conformal gravity [9] and various for- 
mulations of the modified Newtonian dynamics (MOND) 
paradigm [10-13] like, for instance, tensor-vector-scalar 
theory (TeVeS) [14] and quasi-linear MOND [15]. 

Alternative frameworks of the latter kind predict im- 
portant deviations from Newtonian mechanics in regions 
where the total gravitational acceleration is much smaller 
than the value of ao » 10~^° m s~^, i.e. ^jv <C ao- In 
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most parts of the SS, the total gravitational accelera- 
tion exceeds this value by several orders of magnitude. 
As first noted in Ref. [16], however, there exist numer- 
ous points in the SS where the sum of all gravitational 
fields exactly cancels out. These are the saddle points 
(SPs) of the Newtonian potential near the planets which 
are close to, but do not coincide with the Lagrangian LI 
points (where also the centrifugal potential is taken into 
account). These SPs are embedded into ellipsoidal low- 
acceleration regions (bubbles) where the gradient of the 
gravitational potential becomes very small, i.e. compara- 
ble to or smaller than ao, providing an ideal environnic;nt 
for testing some of the predictions of this set of modified 
gravity theories. 

The LISA Pathfinder space mission (LPF) [17] is an 
ESA project currently planned for launch in 2014. LPF 
will test technologies that will be employed in the follow- 
ing LISA (Laser Interferometer Space Antenna) mission, 
which aims at the detection of gravitational waves. The 
LPF will consist of two closely spaced test masses, whose 
distance relative to each other will be measured by means 
of an interferometer. The spacecraft will be shielded from 
any non-gravitational influence such that the masses in- 
side the spacecraft will be in an almost perfect geodesic 
motion. The LPF will therefore be able to measure tidal 
stresses to an unprecedented accuracy. The expected tra- 
jectory of LPF will be a Lissajous orbit around the LI 
point. With minimum corrections to its planned orbit, 
LPF may be redirected to the low-gravity region encas- 
ing the Earth-Sun or Moon-Sun SPs [18]. In previous 
works, it has been claimed that the LPF instruments have 
enough sensitivity to detect the extra tidal stress signal 
predicted by MOND and TeVeS [16, 18]. As such, this 
particular mission is believed to provide a unique occa- 
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sion to put such modified gravity frameworks to a direct 
test. Tire strength of tlic MONDian perturbations that 
the spacecraft will experience when crossing the Earth- 
Sun low-gravity bubble, however, depends strongly on 
the transition speed between the Newtonian and MOND 
regimes, as well as on the intrinsic structure of whichever 
MOND formulation is used. The full space of possible re- 
sults has not yet been explored in detail. In the present 
contribution, we use the recently formulated quasi-linear 
MOND framework (QMOND) [15] to make predictions 
about the expected deviations from Newtonian gravity 
that would occur inside the low-gravity region encasing 
the Earth-Sun and Moon-Sun SPs. Since the positions 
and sizes of the low-acceleration bubbles depend on the 
specific configuration of the SS at a certain instant, we 
also address the problem of determining their positions 
to very high precision. This is a necessary prerequisite to 
any realistic performance of a saddle-fly by experiment. 

The present paper is structured as follows: In Sec. II, 
we introduce the basic concepts regarding QMOND and 
low-acceleration bubbles. In Sec. Ill, we present an accu- 
rate model of the SS Newtonian potential which is used 
to compute the positions and three-dimensional motions 
of the Earth-Sun, Moon-Sun and Saturn-Sun SPs. Fol- 
lowing the lines of previous work [16, 18], we then use our 
Newtonian model to compute the QMOND phantom DM 
distribiition within the low-gravity bubbles surrounding 
the SPs. In Sec. IV, we present and discuss numeri- 
cal solutions to the QMOND potential near the Earth- 
Sun SP. We further derive the extra tidal stress felt by a 
spacecraft near the saddle and readdrcss the question of 
detecting such signals with the LPF. Considering the SP 
between two massive bodies, we also derive semi- analytic 
solutions for the potential which, together with a compar- 
ison to a full numerical treatment, are separately given 
in an appendix. Finally, we conclude in Sec. V. 



II. LOW ACCELERATION BUBBLES IN THE 
SOLAR SYSTEM 

A. Quasi-linear MOND 

QMOND or (QUMOND) [15] is a recently proposed 
quasi-linear formulation of the original MOND paradigm 
[10 12]. Within this framework, the gravitational poten- 
tial $ generated by a mass distribution pbary is a solution 
of Poisson's equation for the modified source density 

A$ = A'kGp = -V • [i^(5Jv/ao)gAr] , (1) 

where ^jv = |gjv| is the Newtonian acceleration and 
the acceleration scale ~ 1.2 x 10~-^° m/s~^ is Mil- 
grom's constant. Here the free fimction i'{y) is related 
to the MOND interpolation function ji{x) according to 
v{y) = l/^{x) and y = xii{x), where x = g/ao is the 
MOND total gravity measured in units of [19]. In 
the deep-MOND regime, i.e. <C ao, we have that 



v{gN /o^o) ~ {gN /cLo)~^^'^ whereas 1 for strong grav- 
itational fields, i.e. gN ^ ao, in which case Eq. (1) 
reduces to the usual Poisson equation. For later pur- 
poses, it is convenient to express the modified QMOND 
density p as the sum of the baryonic contribution pbary 
and a phantom dark matter (PDM) term Ppdm: 

P = Pbary + PPDM- (2) 

As has been discussed in Ref. [15], QMOND is derivable 
from an action and thus satisfies all the usual conserva- 
tion laws. Furthermore, it is important to note that, as 
has been shown in Refs. [20, 21], our Eq. (1) is the gen- 
uine nonrelativistic limit of a particular class of bimetric 
gravity theories [22] . Contrary to TeVeS and similar pro- 
posed relativistic extensions, these constructions do not 
necessarily require a renormalization of the gravitational 
constant, meaning a difference between the bare value of 
G (present in the Lagrangian) or Gc (the effective value 
used in cosmology calculations) and the usual Newtonian 
gravitational constant measured on Earth. In Sec. II B, 
we will further elaborate on this and its implications for 
the present analysis. 

The practical advantage of the QMOND formulation 
with respect to classical MOND is the fact that it re- 
quires only solutions to linear partial differential equa- 
tions, with a further nonlinear algebraic step. In general, 
the QMOND density p can be interpreted as the mod- 
ified source density that yields a MOND-like potential. 
Using simple algebra, we can express Eq. (1) as 

P = (5Ar/«o) Pbary - (47rGao)~^ v'VgN ■ gAf. (3) 

Near a SP (in the absence of baryonic matter) , we have 
Pbary = and the above reduces to 

/5 = ppDM = - (47rGao)~ V'V^jv • gjv, (4) 

where 



is the derivative of the interpolation function. The PDM 

concept is useful to visualize the effects of nonlinear- 
ity in MOND using the standard Newtonian picture: 
For example, it was shown that, at galactic scales, the 
PDM may generate an effective DM disk [23, 24]. Other 
work demonstrated that placing a galaxy into an external 
graviational field leads to an hourglass-shaped distribu- 
tion of negative PDM density [25]. The symmetry axis 
of this hourglass-shaped distribution is parallel to the 
direction of the external field (see, e.g.. Fig. 3 of Ref. 
[25]). Presuming that DM particles always give rise to 
a positive dynamical density, this points toward one of 
the most important differences between (cold) DM and 
MOND at galactic scales. At solar system scales, the sit- 
uation is quite similar since, as we show in Sec. IV, an 
analogous distribution of negative PDM is present near 
the SPs of the total Newtonian potential in the SS. 
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B. QMOND bubbles 

It has been claimed that a direct test of Newton's laws 
in the regime of ultra-low accelerations is feasible with 
our current technological capabilities. More precisely, 
it has been pointed out that a spacecraft, closely ap- 
proaching one of the low-gravity bubbles near the Earth, 
should be able to detect an additional tidal stress signal 
in MOND whose amplitude and shape has been predicted 
in the framework of TeVeS [16, 18]. Naively, one would 
expect these bubbles to be relatively small since MOND 
effects should become important only for accelerations 
smaller or approximately equal to oq. A simple Newto- 
nian calculation shows that the semi-major axis of the 
ellipsoid defined by the condition gjv = aq is w 5 m for 
the Earth-Sun saddle and 1 m for the Moon-Sun sad- 
dle. In TeVeS, however, the dynamics of MOND enters 
through an additional scalar field (f>, i.e. $ « + (f>, 
where $Ar is the usual Newtonian potential (see App. C). 
Thus, even if the total potential $ is still dominated by 
its Newtonian contribution, it is possible to have ^ al- 
ready in the deep-MOND regime. This peculiarity has 
been exploited in Ref [16] to estimate a characteristic 
acceleration scale in TeVeS below which full MONDian 
behaviour of (p should be triggered, OtHg = (47r/fc)^ao, 
where k is the TeVcS scalar parameter. Dividing atiig 
by the tidal stress at a generic saddle point then yields a 
corresponding TeVeS bubble with major axis 
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where T is the Newtonian tidal stress along the symmetry 
axis in the two-body approximation. 
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dsJmjM is the distance of the SP from the lighter 
body (assuming M » m), and d denotes the separation 
between the two bodies. For example, the Newtonian 
tidal stress at the Earth-Sun SP is T « 4 x 10"" s'^. 
Fixing k = 0.03 throughout this work, we have Otrig ~ 
10'''ao and the above formulae allow one to estimate a 
bubble size of 766 km for the Earth-Sun bubble, 280 km 
for the Moon-Sun, and 3.57 x 10^ km for the Saturn-Sun 
bubble. 

As for the framework of QMOND and its rclativistic 
bimetric formulation, we assume that the value of Gc, 
i.e. the effective gravitational constant that is relevant 
to cosmology (e.g. in nucleosynthesis calculations), is 
identical to both the bare G (entering the theory's La- 
grangian) and the value of G measured on Earth (i.e. 
G = Gc = Gbare) [26, 27]. In the context of Milgrom's 
general bimetric MOND gravity [20, 21], such a situation 
explicitly corresponds to setting a + l3 = and (3 = 1. 
This leaves substantial freedom for choosing viable inter- 
polation functions vly) (cf. Ref [28]), and also means 
that our Eq. (1) in Sec. II A is rigorously correct. Despite 



a first investigation of Friedmann-Robertson- Walker so- 
lutions [29], the full cosmological scenario in the frame- 
work of these theories remains to be worked out com- 
pletely, and there are presently no constraints on how 
fast i'{y) may approach the Newtonian limit. In the rest 
of the paper, we therefore take G = Gc = Gbaro, bearing 
in mind that new cosmological constraints on ^{y) might 
arise in the future. An important consequence of this 
assumption is that a characteristic bubble size akin to 
Eq. (6) does not exist in QMOND. Indeed, the bubble 
sizes associated with a detectable anomalous stress sig- 
nal will sensitively depend on the assumed form of v{y), 
a problem we will address in Sec. IV. However, we will 
nevertheless use rtrig given by Eq. (6) to define bubble 
sizes when addressing the problem of their motions and 
evolution in Sec. III. 

A necessary prerequisite for any realistic performance 
of the LPF flyby experiment is the precise knowledge of 
the positions and motions of the low-acceleration bub- 
bles near the Earth. Neglecting extragalactic gravita- 
tional fields (which typically produce accelerations of 
^ext << 10~^° m s~^), the total Newtonian gravitational 
potential at a point within the SS is due to the contri- 
butions of the Sun, the planets (as well as their satel- 
lites) and the gravitational field of the Milky Way. The 
exquisite knowledge of the relative magnitudes of their 
Newtonian potentials allows one a precise determination 
of the SP positions within the SS and the shape of the 
isopotential surfaces nearby. An approximate estimate 
of perturbations induced by a constant external acceler- 
ation field of magnitude Gp on the position of a SP is 
given by (see Ref. [16]) 
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Substituting values for the Earth-Sun SP, we find that 
the Moon induces a mean perturbation on the two-body 
position of the SP of the order of 5r k, 6000 km while the 
presence of Jupiter leads to a perturbation of 5r w 10 km. 
However, the gravitational field generated by the Milky 
Way also plays a role in determining the exact position 
of the SP. Its gravitational acceleration across the SS can 
be estimated assuming that the Sun's orbit aroimd the 
center of the Galaxy is in centrifugal equilibrium. Using 
this hypothesis, we have u^/r « 1.9 x 10~^° m s~^ (as- 
suming V = 220 km s~^, r = 2.5 x 10^^ km), three orders 
of magnitude less than the gravitational pull exerted by 
Jupiter. This would translate into a positional shift of 
the SP of a few cm, and can therefore be neglected. 

From these back-of-the-envelope calculations, we con- 
clude that, within the three-body approximation the po- 
sition of the Earth-Sun SP can be tracked with an ac- 
curacy of «10 km. A plausible impact parameter for 
the LPF spacecraft has been estimated as 6 « 50 km 
[18]. As we will show in the next sections, however, the 
strength of the MOND tidal stress signal decreases sub- 
stantially with the "rapidity" of the transition between 
the Newtonian and MOND regimes, which is determined 
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by the slope of ^{y). As a consequence of this, to test 
a wider range of interpolation functions, the LPF space- 
craft might need to approach the SP more closely than 
previously thought. This motivates our choice to use a 
many-body code to calculate the true motion of the SPs 
to excellent accuracy given by our current knowledge of 
the SS. 



III. MOTION OF THE MOND BUBBLES 

A. The computational framework 

We have used the symbolic calculus package Mathe- 
matica to build a high accuracy model of the SS's New- 
tonian potential and to compute the positions of the 
Earth-Sun, Moon-Sun and Saturn-Sun SPs. In com- 
puting the total Newtonian potential across the SS, our 
model takes into account the planetary potential and 
the contributions of the most massive natural satellites. 
For the positions of the planets, we use INPOP08 high- 
precision ephemeris [40] and adopt the most recent lAU 
estimates for the Gaussian gravitational constant, the 
planet's masses relative to the Sun and their gravitational 
moments. All adopted values and their corresponding 
references arc listed in Table I. Furthermore, note that 
the oblateness of the Sun, Jupiter, Saturn and the Earth 
are taken into account when computing the potential. 
Assuming polar coordinates, the Newtonian potential of 
a planet (centered on the planet with the equatorial plane 
at ^; = 0) is represented by the series 
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where J„ the n-th gravitational moment of the planet, 
r = ^J + + and P„ the Legendre polynomial of 
order n. In our code the series is truncated at n = 2. 
As the distance of the other planets from the SPs is very 
large, for them also the terms of 0(r^^) can be neglected. 
Once all contributions to the Newtonian potential have 
been added together, our code calculates the positions of 
the SPs by numerically solving the equation V(^jv = 
using Newton's tangents method. 



B. The Earth-Sun SP 

In general, the SS SPs of the total Newtonian poten- 
tial follow non-inertial trajectories. An intuitive picture 
of the characteristics of this motion can be obtained by 
starting with the two-body approximation and consider- 
ing the effects of a perturber later. Consider two bodies 
of masses M and m, respectively, with M 3> m. Choose 
a corotating Cartesian coordinate frame centered on the 
orbiting object with the x-axis pointing toward the heav- 
ier body and call the .xy-plane the ecliptic plane. In this 
reference frame, the SP is stationary if the orbit of the 
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FIG. 1. Distance of the Earth-Sun SP from the Earth during 
a year: The gentle yearly variation is due to the eccentricity 
on the Earth's orbit and its ampUtude is (measured on the 
picture) around Argp = 8800 km. The vertical peaks are 
due to the moon's perturbation. Their measured amplitudes 
oscillate between 4700 km < o < 9000 km. 
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FIG. 2. Epicycles described by the Earth-Sun SP around its 
unperturbed position in a year (projected on the plane of the 
ecliptic; here the a;j/-plane). 



lightest body is circular, and its distance from the center 
of the orbit is fixed at rgp. If the orbit is instead ellip- 
tical, then a periodic oscillation of the SP along the line 
joining the two bodies arises. This oscillation has the 
same period as the orbiting body and an amplitude that 
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TABLE I. Table of adopted constants: All values are expressed in SI units. 



Name 


Symbol 


Value 


Uncertainty 


Reference 


Gaussian gravitational constant 


ka 


0.01720209895 


~ 


lAU 1976, [30] 


Astronomical unit 


AU 


149597870700 m 


±3 m 


[30] 


Heliocentric gravitational constant 


GMq 


1.32712442099x10^° m'^s"^ 


±1x10 


[31] 


Geocentric gravitational constant 


GM(s 


3.986004418x10^" m^s"^ 


±8x10® 


[32] 


Mq / MMercury 




6.0236x10® 


±3x10^ 


[33] 


Me/Mvenus 




4. 08523719 xlO'^ 


±8x10"* 


[34] 


Me/MMars 




3.09870359 xlO'^ 


±2x10"^ 


[35] 


Mo/Mjupiter 




1.047348644x10^ 


±1.7x10"® 


[36] 


Mo/Msaturn 




3.4979018x10* 


±1x10"* 


[37] 


Mo/Muranus 




2.290298x10* 


±3x10"^ 


[38] 


M0/MNeptune 




1.941226x10* 


±3x10-2 


[39] 



is given by 

Arsp = Ad^-, (10) 

where Ad is the difference between the apoapsis and pe- 
riapsis of the orbiting body. For the Earth-Snn SP, the 
above equation gives Argp ~ 8666 km. If a perturber is 
added to this two-body system (for example, the Moon 
in the case of the Earth-Sun SP) , then, besides this gen- 
tle yearly variation due to the orbital eccentricity, the 
SP will also describe an epicycle around its unperturbed 
position. The semi-major axis of this epycicle is given by 
Eq. (8). The minor axis of this epicycle is half of its ma- 
jor axis and it is parallel to the line joining the two major 
bodies. The period of this epicycle is given by the period 
of the perturber (if it orbits around the lighter body). 
For the Earth-Sim SP, the perturber is the Moon. The 
mean semi-major axis of the Earth-Smi bubble epicycle 
due to the Moon's perturbation is a ~ 6000 km and its 
period equals the synodic period of the Moon. Since the 
Moon orbits on a plane inclined with respect to the eclip- 
tic plane, the SP has also a component of motion in the 
direction perpendicular to the ecliptic plane. Therefore, 
the yearly motion of the Earth-Sun SP is a complicated 
three-dimensional trajectory, a composition of the differ- 
ent dynamical effects mentioned above. Using our code, 
we have produced several plots to illustrate these com- 
plicated motions. In Fig. 1, we show the distance of the 
Earth-Sun SP from the Earth during a year. 

The gentle yearly variation of the distance due to the 
Earth's eccentric orbit is roughly Argp = 8800 km, very 
close to the value calculated with the help of Eq. (10). 
The vertical spikes caused by the Moon's presence have 
amplitudes varying between 4700 km< a <9000 km, with 
a mean value of around 6850 km. This is compatible 
with the amplitude calculated from Eq. (8). In Fig. 2, 
we show the epicycles (projected onto the ecliptic plane) 
described by the Earth-Sun SP around its unperturbed 
(two-body) position. In this plot, the Earth lies on the 
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FIG. 3. Distance of the Earth-Sun SP from the ecliptic plane: 
This component of the motion is due to the inclination of the 

Moon's orbital plane with respect to the ecliptic. The peak- 
to-valley amplitude of this yearly oscillation is around 5000 
km. 

right-hand side. At new moon, the SP shifts toward the 
Earth, whereas it shifts toward the Sun at full moon. 
The epicycle's amplitude is higher if new moon occurs 
at perigee and lower if at apogee. The yearly motion of 
the Earth-Sun SP in the direction perpendicular to the 
ecliptic plane is illustrated in Fig. 3. 



C. The Moon-Sun SP 

Another SP close enough to the Earth for direct testing 

is the Moon-Sun saddle [18], sometimes also referred to 
as the "Earth- Moon" SP [16]. However, since the gravi- 
tational pull of the Sun on the Moon is higher than that 
exerted by the Earth, it is better to look at it as a Aloon- 
Sun SP heavily perturbed by the Earth [18]. The magni- 
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FIG. 4. Distance of the Moon-Sun SP from the Moon during 
a sidereal month: The asymmetry of the bell-hke curve is due 
to the inclination of the Moon's orbit. 



tude of the perturbation induced by the Earth on the po- 
sition of the bubble is Arsp w 13000 km, where wc have 
again used Eq. (8). In Fig. 4, we plot the Moon-Sun 
SP distance from the Moon for a whole sidereal period. 
The epicycle's of the SP around its unperturbed position 
is shown in Fig. 5. The measured variation in distance 
during a sidereal month is Argp « 12550 km which agrees 
well with the calculated value. The SP reaches its mini- 
mum distance from the Moon at full moon (when it lies 
between the Earth and the Moon) , and its maximum dis- 
tance during new moon (when it lies between the Moon 
and the Sun). The inclination of the Moon's orbit with 
respect to the ecliptic changes the Moon-Sun saddle's 
apoapsis. This gentle variation can be seen in Fig. 6, 
where we show the SP distance from the Moon's cen- 
ter for a whole year. The peak-to-valley variation in the 
apoapsis has an amplitude of approximately 2000 km. 

A question raised in Ref. [16] was whether it may 
happen (considering the large magnitudes of perturba- 
tions) that the Moon-Sun SP crosses the Moon's surface 
during its motion. Prom our code, we find a minimum 
distance of the saddle from the Moon's center of 22512 
km. Since the Moon's radius is 1738 km, we conclude 
that the SP never crosses the Moon's surface. In Fig. 
7, we show the motion of the SP in the direction per- 
pendicular to the ecliptic (due to the inclination of the 
Moon's orbit). For the Moon-Sun saddle, the amplitude 
of this motion component is lower (3000 km compared to 
5000 km) than in the case of the Earth-Sun saddle. This 
can be explained as follows: Whereas the Moon (with his 
inclined orbit) tends to push the Earth-Sun saddle away 
from the ecliptic plane, the combined pull of the Sun and 
the Earth tends to push the Moon-Sun SP toward the 
ecliptic plane. 
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FIG. 5. Epycicles described by the Moon-Sun SP around its 
unperturbed position during a year: Here the Moon lies on 
the left-hand side. 



D. The Saturn-Sun SP 

The Saturn-Sun low-acceleration bubble is one of the 
largest in the SS. The semi-inajor axis of the MOND 
bubble is 3.57 x 10*^ km k, 0.02 AU. The minimum dis- 
tance between Jupiter and Saturn is around 5.3 x 10^ km. 
When this minimum is reached, Jupiter exerts a gravita- 
tional pull of flp w 4.4 X 10~''m s~^ on the Saturn-Sun 
SP. The Newtonian two-body tidal stress at the Saturn- 
Sun SP calculated from Eq. (7) is T w 5.9 x IQ-i^ 3-2, 
Using Eq. (8), we may therefore calculate the peak-to- 
valley amplitude of this perturbation which turns out to 
be around 1.5 x 10^ km. This is much lower than the vari- 
ation in distance due to the orbital eccentricity which is 
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FIG. 6. Distance of the Moon-Sun SP from the Moon over a 
year. 




FIG. 7. Distance of the Moon-Sun SP from the plane of the 
ecliptic: This component of the motion is due to the incli- 
nation of the Moon's orbital plane with respect to the eclip- 
tic. The peak-to-valley amplitude of this yearly oscillation is 
around 3000 km. 



approximately 2.54 x 10^ km. In Fig. 8, we show the 
distance of the SP from Saturn during a Saturnian year. 

Considering the size of the Saturn-Sun MOND bub- 
ble, there is a considerable chance that minor bodies or 
some of Saturn's outer satellites periodically intersect it. 
Indeed, we have found that a subsample of the Norse 
group [41, 42], consisting of Surtur, Ymir, Fenrir and 
Loge, periodically go through the Saturn-Sun MOND 
bubble. These small satellites arc characterized by retro- 
grade motion. The semi-major axes of their orbits vary 
from 2.24x lO"^ km (Fenrir) to 2.30x 10^ km (Loge). Their 
orbits are periodically intersected by the MOND bubble 
whose distance from Saturn varies between 2.24 x lO'' 
km to 2.5 X 10'' km (see Fig. 8). Furthermore, one finds 
that the satellites' inclination with respect to the eclip- 
tic is low and varies from 3° (Surtur) to 16° (Fenrir). 
Since the MOND bubble stays always close to the ecliptic 
plane, close encounters with these objects are generally 
favored. An example of such an encounter between Fen- 
rir and the Saturn-Sun bubble is shown in Fig. 9. For all 
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FIG. 8. Distance of the Saturn-Sun saddle point from saturn 
during a Saturn's year: The amplitude of the mainly due to 
the eccentricity of Saturn's orbit is around 0.017 AU. 




2000 4000 6000 8000 II 

Time [days] 

FIG. 9. Distance of the Saturn's moon Fenrir from the Saturn- 
Sun SP during a Saturnian year: The horizontal dashed line 
shows the size of the MOND bubble according to Eq. (6). 
Note that 6212 days after the 2012 Jan 1, the Moon gets very 
close to the center of the bubble at a distance of only rmin 
0.437rtrig w 1.43 x 10^ km. 



of these satellites, we list their dates of encounter with 
the Saturn-Sun bubble, their minimum distances from 
the SP as well as their crossing times in Table II. 

So far, the found subset of Saturn's satellites comprises 
the only known objects in the nearby SS that periodically 
enter the MOND regime. Even though the MONDian in- 
fluence on the satellites' trajectories is likely to be very 
small, it cannot be ruled out that the induced perturba- 
tions at each encounter may accumulate to a potentially 
observable effect. Regarding a detailed analysis of the 
MONDian impact on the path of these satellites, we re- 
fer to a forthcoming paper. 



E. Error estimates 

The most important sources of error on the computed 
positions of the SPs are the uncertainties on the plan- 
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TABLE II. Norse's satellites minimum distances to the 
Saturn-Sun MOND bubble. 



Satellite 




dmin 1 '^trig 


Date; 


Crossing t 




(10** km) 




(dd/mm/yyyy) 


(days) 


Surtur 


2.00 


0.48 


24/06/2015 


64 


Ymir 


0.65 


0.16 


09/02/2013 


68 


Fenrir 


1.43 


0.43 


03/01/2029 


49 


Loge 


0.60 


0.17 


16/03/2026 


58 




400 600 
Iterations 



1000 



FIG. 10. Monte-Carlo simulation output: This plot shows the 

distance of the perturbed Earth-Sun SP from its unperturbed 
position. The mean value of this data series is a good estima- 
tor of our final uncertainty on the position of the Earth-Sun 
SP. 



ets' positions and masses. We have used a Monte-Carlo 
method to estimate uncertainties in the positions of the 
SPs. The masses of the planets are perturbed using ran- 
dom values drawn from normal distributions with means 
at the quoted parameter values, and standard devia- 
tions matching those listed in Table I. For the error on 
the planets' position we use the uncertainties in the IN- 
POP08 ephemeris [40] which we assume to be roughly « 

3 m for the Earth's position, a few cm for the Moon's po- 
sition, « 1000 km for the giant planets and « 10 km for 
the inner planets. With this method, we obtain an abso- 
lute error on the position of the Earth-Sun SP of Ae-s ~ 

4 m (see Fig. 10). This quantity is close to the error on 
the two-body position of the Earth-Sun SP Ag_^°'^'*'^ w 

5 m which we obtained by propagatin g the e rrors quoted 
in Table I to the quantity rgp k, d^Jm/M representing 
the approximate distance of the Earth-Sun SP from the 
Earth-Moon system barycenter. Our results tell us that 
the major sources of error on the position of the Earth- 
Sun SP are the uncertainties on the GMq/GMq ratio 
and on the AU. The contribution of the other planets 
to the final error sum up to an amplitude of about 1 m. 
For the Saturn-Sun SP we obtain, in the same way, an 
uncertainty of As at- S ~ 2 km. 



IV. QMOND POTENTIAL INSIDE THE 
BUBBLE 

A. Phcintom dark matter density necir the SP 

As can be seen from Eq. (4), the PDM density inside 

a MOND bubble is entirely determined by the Newto- 
nian potential and the choice of a specific interpolation 
function. Using our model of the SS gravitational po- 
tential, we have calculated the PDM density distribution 
around the Earth-Sun, Moon-Sun and Saturn-Sun SPs 
for different choices of v{y). In particular, we have used 
the QMOND analog of the so-called simple interpolation 
function, ii{x) = x/{x + 1), which takes the form 



(11) 



This functional form has been shown to provide a good 

description of extragalactic phenomenology and allows 
one, for instance, to fit the observed rotation curves of a 
wide range of spiral galaxies very well [43]. In addition, 
we have considered the interpolation function 



1/4 



(12) 



which is the QMOND analog of the function used in Refs. 
[16, 18] and exhibits a behavior close to that of the func- 
tional form originally proposed by Bekenstein in the con- 
text of TeVeS [14]. Finally, we have also chosen an inter- 
polation function similar to Eq. (11), but giving rise to 
a "faster" transition between Newtonian dynamics and 
MOND: 



16 ) 



(13) 



This choice is motivated by recent claims that interpo- 
lation functions like Eq. (11), which still result in a 
relatively "slow" transition between the two dynamical 
regimes, could give rise to an anomalously high secular 
precession of the outer planet's perihelion, incompatible 
with published residuals [44, 45]. In Fig. 11, we plot the 
three interpolation functions as a function of the gravi- 
tational field strength expressed in units of uq. As can 
be seen, the function defined by Eq. (12) approaches the 
value I'd/) = 1 (corresponding to the Newtonian regime) 
much more slowly than the other ones. Considering their 
performance on galactic scales, we briefiy discuss all three 
interpolating functions in App. B where it is also demon- 
strated that the behavior of Eq. (13) at these scales is 
very similar to that of the simple z^-function Eq. (11). 

As remarked earlier, the PDM density depends on z/', 
which basically tells us that the PDM density within the 
MOND bubble is sensitive to the "transition" speed be- 
tween the dynamical regimes. In fact, choosing an inter- 
polation function with a fast enough transition, it is pos- 
sible to suppress non-Newtonian effects for accelerations 
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>' = wA'o 

FIG. 11. Interpolation functions as a function ofy — gjv/ao: 
Solid, dotted and dashed lines correspond to Eq. (12), Eq. 
(ll)andEq. (13), respectively. Vertical dashed lines show the 
characteristic scales J/ = 1 and y = (47r/fc)^ for k — 0.03. The 
Bekenstein-like interpolation function Eq. (12) gives rise to a 
much slower transition between the two dynamical regimes. 
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FIG. 12. Logarithmic plot of the rescaled PDM density's 
absolute value 47rG/9pDM at a cut of j/ = 50 km near the Earth- 
Sun SP: The line styles are defined as in Fig. 11. Note that 
the PDM density decreases by several orders of magnitude for 
the (/-functions with increased slope. 



above ag. The PDM density near the SP is proportional 
the trace of the stress tensor, 47rG'ppDM = tr(5'ij), where 
Sij = didj^, and provides a good order-of- magnitude 
estimator for the maximum MOND tidal stress signal 
within a certain region - if one excludes values at special 
points such as the SPs themselves. 

The above dependency on the slope of the interpolation 
function can be seen in Fig. 12 where we plot the PDM 
density at cuts of y = 50 km near the Earth-Sun SP. It is 
evident that the PDM density changes by several orders 
of magnitude between the different choices of v{y). We 



FIG. 13. PDM isodensity contours inside in a box of 5 x 10^ 
km centered on the Earth-Sun SP along the z = Q plane and 
assuming v given by Eq. (12): The positions of the Earth and 
the Moon-Sun SP are marked with a large and a small dot, re- 
spectively. The arrow points towards the Sun and the dashed 
line indicates the ppdm = contour. The PDM density is 
negative within a conic region (deformed by the Moon's pres- 
ence) , with the symmetry axis perpendicular to the Earth-Sun 
direction, and positive elsewhere. The Earth is embedded into 
a halo of positive PDM. Moreover, the PDM density gradient 
is bigger near the Moon-Sun SP. 



therefore expect that at the same distance (or passing 
trajectory) from the SP, the corresponding tidal stress 
signal will be much smaller for an interpolation function 
given by Eq. (11) than for Eq. (12). We will confirm 
this preliminary prediction by using numerical solutions 
of the QMOND potential near the Earth-Sun SP below. 
In Fig. 13, we illustrate the PDM density distribution 
calculated with the help of Eq. (12) and rescaled by a 
factor of 4 X IO^^ttG within a box of 10^ km edge length 
centered on the Earth-Sun SP. As expected, the PDM 
density distribution around the SP exhibits cylindrical 
symmetry, with the symmetry axis along the line joining 
the planet with the Sun. In the case of the Earth-Sun 
SP, this symmetry is perturbed by the presence of the 
Moon. 



B. Multi-grid FFT solver 

Since the QMOND PDM decreases quickly toward zero 
when moving away from the SP, the boundary value 
problem given in Eq.(l) may be solved by means of effi- 
cient numerical methods based on the Fast Fourier Trans- 
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form (FFT). To achieve this, we have implemented a 
modified version of the Poisson solver described in Ref. 
[46] which now incorporates a multi-grid scheme to ob- 
tain a higher resolution near the SP. 

The basic idea is the following: Beginning with a large 
physical box and a coarse sampling of the underlying den- 
sity distribution on an equidistant grid, we solve for the 
potential $ assuming zero boundary conditions - since 
we arc only interested in the purely non-Newtonian con- 
tribution = $ — $jv- Keeping the same number 
of grid points per dimension N, we then move to a grid 
of smaller physical size such that the previous solution 
remains sufficiently accurate at the new boundaries. A 
first correction to the potential is then found by consider- 
ing the solution for the correction density field, which is 
obtained as the difference between the true PDM density 
and the one calculated by interpolating the density field 
on the coarse grid. For this step, we again assume zero 
boundary conditions and use a cubic spline for the in- 
terpolation. Finally, the above procedure is successively 
repeated until the desired resolution in the central parts 
is reached, and adding the individual results yields the 
full solution with the appropriate boundary conditions. 
Apart from the usual sanity checks (see, e.g., Ref. [46]), 
the method was applied to several input densities and 
the corresponding results were directly compared to those 
obtained with the original single-grid solver. For our pur- 
poses, we have found that the relative differences between 
these results are basically negligible, typically on the or- 
der of < 10~^. In App. A, we also present semi-analytic 
solutions (both for the quasi-Newtonian and the deep- 
MOND regime) which are used for comparison to the 
full numerical results. 

In the case of the Earth-Sun saddle, for instance, we 
set N = 448 and use three physical cubic boxes with sides 
of 20000, 5000 and 1250 km, respectively. This yields a 
maximum resolution of approximately Axmax = 2.8 km 
in the saddle region. Note that all grids are chosen such 
that the SP is located at (0, Aa;max/2, 0) to avoid any 
numerical difficulties arising from the divergence of the 
PDM density at the exact SP. Assuming these parame- 
ters, the numerical solution is found to be reliable for dis- 
tances > 20 km from the SP (again see App. A). Compu- 
tationally, our spectral multi-grid method performs quite 
well: Considering the above setup, the full solution and 
derived quantities such as the stress tensor are calculated 
on a timescale of less than half an hour on an average ma- 
chine using a single processor. 



C. Numerical results 

Using the above multi-grid FFT code, we have numeri- 
cally solved Eq.(l) near the Earth-Sun saddle point. The 
input PDM densities have been calculated for the two- 
body (Earth-Sun) and the three-body (Earth-Sun-Moon) 

approximations. Hereafter, we adopt a reference frame 
centered on the SP, with the x-axis pointing along the 
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FIG. 14. Two-body QMOND tidal stress signal for « = 
along the lines y w 25 (solid line), y w 100 (dashed line) and 

y « 400 km (dotted line) inside the Earth-Sun MOND bubble: 
For the calculation, we have assumed the interpolation func- 
tion Eq. (12). At y ~ 25 (?/ ~ 400) km, the positive (negative) 
peak is at Syy « 7.5 x lO-^s'^ {Syy « -0.2 x IQ-^s"^). 



Earth-Sun direction (increasing toward the Sun). Note 

that within this reference frame, which is identical to the 
one used in Ref. [18], the LPF instrument will only be 
sensitive to the Syy component of the stress tensor be- 
cause of certain design characteristics [47]. The results 
of our simulations for the Earth-Sun bubble are shown 
in Figs. 14, 15 and 16: Using the Bekenstein-like inter- 
polation function given by Eq. (12), Fig. 14 shows the 
two-body QMOND tidal stress signal for 2; = at cuts 
of around y = 25, 100 and 400 (±2.8 due to the fixed 
sampling on equidistant grids) km away from the Earth- 
Sun SP. Our result fairly reproduces that of Ref. [18] (cf. 
their Fig. 6). Compared to the calcidation for TcVeS, the 
main differences are a generally larger tidal stress signal 
at all cuts (by about a factor of 1.7) and a deeper central 
wiggle. Note that this agrees well with the conclusions of 
Ref. [28] about "Type 11" MOND-like theories and the 
softening effect of the TeVeS curl field. 

We have also computed the tidal stress signal using the 

simple interpolation function Eq. (11). As expected from 
the earlier analysis of the PDM density, the tidal stress 
signal decreases by approximately two orders of magni- 
tude in this case. From Fig. 15, one finds that the peak 
value of the tidal stress at ^ = and y « 25 km is ap- 
proximately 8 X 10~^®s~^. Including the Moon into the 
calculation (see Fig. 16) has only a small impact on the 
estimated stress signal. We have also verified that the 
tidal stress signal becomes even lower {Syy w lO^-'^^s"^ 
at 2; = and y « 25 km) if the "fast-transition" interpo- 
lation function Eq. (13) is used. In the next section, we 
shall discuss the predicted additional stress signal in the 
light of the expected LPF sensitivity. 



11 




400 



FIG. 15. Same as Fig. 14, but now assuming the simple 
interpolation function Eq. (11): At j/ w 25 (y ^ 400) km, 
the positive (negative) peak is at Syy w 8 x 10~^^s~^ {Syy w 
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FIG. 16. Thrcc-body QMOND tidal stress signal within the 
Earth-Sun bubble for 2 = along the line j/ ~ 25 km for dif- 
ferent lunar phases, computed using the simple interpolation 
function Eq. (11): The solid (dotted) line corresponds to the 
tidal stress signal at approximately full (new) moon. Note 
that the presence of the Moon has only a slight impact on the 
signal. 



D. LPF sensitivity and the QMOND signal 

Choosing the analog of the Bekenstein-like z/-function 
Eq. (12) in the framework of TeVeS, it was concluded 
that the instrument onboard the LPF spacecraft should 
have enough sensitivity to detect the additional MON- 
Dian tidal stress signal [16, 18, 28]. At a distance 
from the SP similar to the predicted impact parame- 
ter, this signal was estimated to be on the order of 



tions, however, the expected signal strength is very sen- 
sitive to the specific form of viy), and one might question 
whether such a conclusion remains valid for other choices. 
In particular, we will be interested in the situation for in- 
terpolation functions like Eq. (11) which provide a good 
description of the phenomenology on cixtragalactic scales. 

Apart from the instrumental characteristics, the sen- 
sitivity of the LPF probe to tidal stresses depends on 
the frequency of the MOND signal [48]. When the LPF 
spacecraft crosses the Earth-Sun bubble, the frequency 
scale of the MOND signal will be given by / = w/Zfe, 
where v is the spacecraft's velocity and denotes the 
characteristic length scale of variations in the MOND 
tidal stress. It can be seen from Figs. 14 and 15 that 
the QMOND tidal stress signal typically varies over a 
length of lii ~ 100 km for both choices of the interpola- 
tion function. Note that this scale length is the same 
as the one adopted in Rcf. [28]. Given a model for 
the power spectral density (PSD) Sh{f) of the noise, 
the expected signal-to-noisc ratio (S/N) may be roughly 
estimated from the signal's characteristic amplitude he, 
S/N ~ hc/y/fSh [49]. To obtain a sensitivity threshold 
for the LPF, we set S/N ^ 1. Considering v = 1.5 
km s~-^ [50] and a constant noise PSD with y/Sh ~ 
1.5 X 10-"s-V\/Hz [28], this gives / « 10"^ Hz and 
yields a threshold on the order of 10~^*s~^. This should 
be sufficient to detect the QMOND tidal stress signal for 
interpolation functions like Eq. (12) whose amplitude 
lies between 2 x 10"" and 8 x 10"" s'^ at j/ = 50 km 
(see Fig. 14). For the simple function Eq. (11), however, 
the situation is different: Since the tidal stress signal de- 
creases by approximately two orders of magnitude, we 
expect that the signal will not be detectable at the given 
impact parameter. Therefore, to have a chance of prob- 
ing these models, the LPF spacecraft might need to pass 
the Earth-Sun SP much closer than previously thought, 
but it is also possible that the modified signal remains 
completely undetectable. 

To get more insight on the above, we follow the lines 
of Rcf. [28] and consider the signal's amplitude spectral 
density (ASD) which is given by the square root of the 
PSD: 



P{f) 



2 
f 



/ h{t) exp{-2Trift)dt 



-T/2 



(14) 
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As we have seen in the previous see- 



where / is the frequency, h{t) = Syy{vt, b, 0) corresponds 
to the signal, and b is the impact parameter, i.e. the 
minimum distance from the SP which is approached at 
t = 0. Fixing v = 1.5 km [18, 28], we furthermore 
have the integration period T which we choose as T w 
1.3 X 10**s, in accordance with our numerical setup (see 
Sec. IV). To compare the signal ASD to that of the noise, 
we use the simplified LPF noise model discussed in Rcf. 
[47] which assumes a constant baseline within 1 < / < 10 
mHz, with an ADS around 1.5 x 10~^^s~^/\/Hz. Beyond 
this range, we assume that the sensitivity degrades with 
1// and for lower and higher frequencies, respectively. 
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FIG. 17. Top panel: The ASD of the MONDian tidal stress 
signal for v(\j) given by Eq. (12) compared to the simple 
noise model described in the text; assuming v = 1.5 km s~^, 
wc illustrate results for 6 ~ 45 (solid line), 25 (dashed lino), 12 
(dotted line) and 1 km (dot-dashed line). Bottom panel: Same 
as the top panel, but now assuming the simple interpolation 
function Eq. (11). 



Note that although this simple noise model will suffice for 
obtaining realistic estimates, improved LPF noise models 
have recently been presented in Ref. [28] (see their Fig. 
6). 

For the interpolation function Eq. (12), the top panel 
of Fig. 17 illustrates the resulting ASDs of the signal, 
assuming different impact parameters between approxi- 
mately 45 and 1 km. Since the numerically calculated 
stress becomes unreliable for distances < 20 km from the 
SP (see App. A), our results for small impact parame- 
ters underestimate the signal at the high-frequency end, 
i.e. for / > 7 X 10~^ Hz. However, this corresponds 
to a frequency domain where the noise is dominant, and 
thus we expect this resolution effect to have no signif- 
icant impact on our conclusions. As can be seen from 
the figure, SjN is typically on the order of 10 per \J Hz 
over a broad frequency range. Applying the techniques 
of noise matched filtering (which are frequently used in 
the search for gravitational waves) to the present prob- 
lem, it can further be shown [28, 49] that the maximal 
SjN , realized by an optimal filtering template, is given 
by 
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Here Sh{f) again denotes the PSD of the noise and h{f) 
is the Fourier transform of the signal. Using Eq. (15), 



the above scenario for the Bckenstein-like interpolation 
bmction Eq. (12) translates into integrated S/N ratios 
of around 50 — 105, which is in good agreement with the 
TcVcS results of Ref. [28]. 

In the bottom panel of Fig. 17, we show the corre- 
sponding signal ASDs for the simple interpolation func- 
tion Eq. (11). Almost throughout the entire frequency 
range, the calculated signals now appear by at least an 
order of magnitude below the noise level. For an impact 
parameter of 6 « 12 km and larger, this gives rise to 
an integrated S/N of around 1 or less. At 6 ~ 45 km, 
for instance, we have S/N « 0.5. Only in the case of 
6 PS 1 km, we find that the maximal S/N crosses unity, 
taking a value of 2. To test models based on interpola- 
tion functions similar to Eq. (11), the LPF spacecraft 
would therefore be required to pass the Earth-Sun SP at 
a distance on the order of a 1 km or less. Note, however, 
that since the expected integrated S/N ratios are rather 
close to unity, these estimates are quite sensitive to the 
details of the true underlying noise model which is cur- 
rently unknown. Depending on the actual noise model, 
this could, in principal, give rise to exposing the modified 
signal to a wider radial range or making it virtually unde- 
tectable. Also, considering potentially necessary impact 
parameters b < 1 km, practical issues related to posi- 
tional accuracies as well as the spacecraft's self-gravity 
might be of concern [28]. Taking all of this into account, 
it seems safe to conclude that the simple interpolation 
function Eq. (11) has a good chance of surviving a LPF 
null result. If the "fast-transition" interpolation func- 
tion Eq. (13) is adopted, then the S/N ratios are even 
further suppressed, and the modified tidal stress signal 
will basically become invisible to the LPF. Finally, note 
that wc have also considered the framework of TeVeS and 
addressed the question whether there exist viable inter- 
polation functions which give rise to similar conclusions. 
For this purpose, we have used our previous choices of 
^{y) to formally construct the corresponding functions /t 
for the TeVeS scalar field. A detailed description of this 
procedure and the implications of our results are sepa- 
rately discussed in App. C. 



V. CONCLUSIONS 

Low- acceleration regions encasing the SPs of the grav- 
itational potential in the SS, may exhibit significant de- 
viations from Newton's laws within the framework of 
MOND or closely related theories. As the LPF space- 
craft may visit the Earth-Sun low-acceleration bubble in 
the near future, this could provide a unique occasion to 
put MOND-like modifications of gravity to a direct test. 
While it is known that the quoted LPF sensitivity is good 
enough to detect the non-Newtonian tidal stress signal in 
certain cases, it is still very important to explore the full 
space of theoretical models and possible outcomes of the 
planned fly by experiment. Also, a necessary prerequi- 
site for any realistic performance of this and future ex- 
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periments, is the precise knowledge of the positions and 
motions of such bubbles near the Earth. 

Using a high accuracy model of the SS Newtonian po- 
tential, we have determined the positions and tracked the 
motions of the Earth-Sun, the Moon-Sun and Saturn-Sun 
SPs with an accuracy of Ae-s ~ 4 m and Asat-S ~ 2 
km for the Earth-Sun SP and Saturn's SP respectively. 
During the analysis of the SP's motion we have further 
found that a subset of Saturn's satellites, consisting of 
Surtur, Ymir, Fenrir and Loge, all regularly pass through 
the Saturn-Sun MOND bubble. This provides the first 
example of SS objects periodically entering the MOND 
regime. An accurate study of the motions of these satel- 
lites, may in principle be used as an independent test of 
the MOND paradigms, but clearly more investigations 
are needed. We plan to perform a quantitative analysis 
on how the MONDian perturbations impact the satel- 
lites' trajectories in a forthcoming paper. 

Adopting the framework of QMOND, we have calcu- 
lated the additional tidal stress signal due to the PDM 
density near the Earth-Sun SP. Using a Bekenstein-like 
interpolation function and assuming an impact parame- 
ter 6 « 50 km, we obtain an approximately by a factor of 
1.7 higher tidal stress signal than estimated in the context 
of TeVeS, which is in accordance with previous considera- 
tions [28]. Choosing the QMOND analog of the so-called 
simple interpolation function, which provides a good de- 
scription of extragalactic phenomenology, we encounter 
a different situation. In this case, the signal turns out to 
be significantly smaller and, taking the sensitivity of the 
LPF probe into account, we typically obtain maximal (in- 
tegrated) S/N ratios on the order of 1 or less. Although 
moving to very small impact parameters on the order of 
^ 1 km might, in principle, help to detect the anomalous 
stress signal, it seems safe to conclude that interpolation 
functions similar to the simple one have a good chance of 
surviving a LPF null result. For interpolation functions 
with a "faster" transition from Newtonian dynamics to 
MONDian behavior, the S/N ratios are even hirthcr sup- 
pressed, and the modified tidal stress signal will basically 
become invisible to the LPF. 

Finally, we want to point out that our conclusions do 
not necessarily apply to other proposed MOND frame- 
works. Considering the original form of TeVeS, for in- 
stance, a transition corresponding to Eq. (13) yields an 
unacceptable interpolation function jl for the scalar field, 
which is briefly discussed in App. C. Moreover, further 
constraints on the rapidity of i^{y) in QMOND might 
arise once the full cosmological scenario of its associated 
relativistic bimetric theory has been worked out. 
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Appendix A: Semi-analytic solutions near the SP 

In the following, we shall consider the configuration 
of two bodies at distance d with masses M and m, re- 
spectively, where we assume that the line connecting the 
bodies coincides with the z-axis. Shifting the coordinate 
origin to the SP and introducing spherical polar coordi- 
nates with z = rcosO, the Newtonian potential in the 
vicinity of the saddle approximately takes the form 



= Cr^ (l - 3cos^ e) + const, 



(Al) 



where 
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If M m (as in case of the Earth-Sun SP, for instance), 
Eq. (A2) may be well approximated by 



(A3) 



Near the SP where A$Ar = 0, the QMOND contribution 
to the total gravitational potential 

$ = $JV + $QM, (A4) 

is associated with the phantom source distribution 

Pqm = A$QM = V • [i/(t/)V$Ar] . (A5) 

Although Eq. (A5) does generally not exhibit analytic 
solutions, it is possible to find approximate expressions 
for certain cases, which we shall discuss in the following 
sections. 

1. QuEisi-Newtonicin domain 

If the resulting dynamics is sufficiently close to the 
Newtonian limit, i.e. y = gN /clq ^ 1, we may expand 
the interpolating function y{y) in powers of y, 



v{y) « 1 + ^ + ^ + . . . . 

V y 



(A6) 



Substituting the above into Eq. (A5) and remembering 
that A$ AT = near the saddle, we thus obtain 
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The linearity of Eq. ( A7) allows one to solve the equation 
order by order. Using Eq. (Al), one finds that the nth- 
order contribution to the source term can be written as 



-(") 
Pqm 



l-9cos2 



r-" (l + 3cos2^)("+')/2' 
where we have defined 



(2C) 



n-l • 



n e N, (A8) 



(A9) 



Let us begin with the first-order contribution. As we 
are only interested in the saddle region, the approximate 
Newtonian potential in Eq. (Al) suggests an ansatz of 
the form 



^om(^^) = Tir-'^^r^W + const, 



(AlO) 



where /3i and 71 are assumed as constant. Inserting the 
above expression into Eq. (A7) then fixes /3i = 1 and 
7i = Al, leaving us with a single ordinary differential 
equation for T'^^(^), 



1 - 9 cos2 B 



(All) 



(l + 3cos2 efl'^' 



For a physically relevant solution, we require that T^^^(0) 
is smooth and periodic on the interval [0,7r] (see, e.g., 
[16]). In this case, Eq. (All) may numerically be solved 
in terms of a Fourier scries, yielding 



r^^^ (6*) « 0.06310 - 0.23719 cos 261 

+ 0.02659 cos 46* - 0.00498 cos f 



(A12) 



The above series converges rather quickly and may al- 
ready be cut off after the first eight terms to achieve a 
relative accuracy of < 10"'*. 

Next, we shall consider the second-order correction to 
the potential $qm- It turns out that one cannot adopt 
the same ansatz as in Eq. (AlO) because it is impossible 
to meet the smoothness condition of the angular part 
in this case. As we will see shortly, this can easily be 
remedied by considering an additional term proportional 
to log(r/rs) in the potential. Thus we start from 



$(2) 



+ const. 



(A13) 



Again substituting this expression into Eq. ( A7) , we find 
that /32 = and 72 = A2. This requires the constant B to 
be dimensionless and gives rise to the following equation 
for r(2)(e): 



cos^ d 



l-9cos2 6i 
(l + 3cos2 9f 



B. 



(A14) 



Since only derivatives of T'-^^O) appear in Eq. (A14), 
one may integrate once and finally obtain 



dO 



7-(2)(^) 



sm( 



2cos6» 

x/3 



BcosO + D 



arctan 



(V3cos6») 



(A15) 



where D is an integration constant. Applying the peri- 
odicity and smoothness conditions to the above, a brief 
calculation leads to I? = and 



B 



V3 



(A16) 



Once again, an approximate solution of T^^-* (9) may be 
found by using a Fourier series approach, and the result- 
ing scries, 



T^^\0) « -0.10460cos26' + 0.02177 cos 40 
- 0.00524 cos 69 + 0.00136 cos + 



(A17) 



quickly converges toward the real solution. In principle, 
one could apply the general procedure outlined in this 
section also to higher orders, but for the purpose of this 
paper, we shall not take these into account. Finally, note 
that the semi-analytic solutions found in this section au- 
tomatically satisfy the desired N(nvtonian boundary con- 
ditions because of Pi, 1^2 < 2, meaning that the total 
gravitational potential is entirely dominated by the ex- 
pression in Eq. (Al) when moving to radii much larger 
than the true size of the bubble, i.e. r/ro » 1, where 
'^o = ''trig for the interpolation function defined in Eq. 
(12). 



2. Deep-MOND domain 

In the limit of very small gravitational fields, we have 
y — )■ and the QMOND interpolating function may be 
approximated as 



1 



(AI8) 



which is independent of the particularly chosen form of 
y{y). Using Eq. (A18) in Eq. (A5), the corresponding 
phantom source term can be expressed as 



PQM = - - 



A l-9cos2 6l 



^''(l-F3cos2 9) 



5/4' 



where 



A=^^/2C^o- 



(A19) 



(A20) 



Similar to our approach in the quasi-Newtonian situa- 
tion, we make the ansatz 



^QMir, 9) = -/r^T{9) + const 



(A21) 
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which is again substituted into Eq. (A5). This time we 
obtain /3 = 3/2 and 7 = A, with the resulting differential 
equation for T{d) given by 



15^,^, cos 6* d <P 
4^(^)+sin?d.^(^) + ^^(^) 
1 - 9 cos2 B 



(A22) 



(1 + 3cos2 



VV4- 



Analogously to the first-order contribution in the quasi- 
Newtonian domain, an approximate solution is deter- 
mined by expanding I'iff) into a Fourier series on the 
interval [0,7r]. A straightforward calculation then yields 



T[B) ~ -0.06750 - 0.53716 cos 26* 

-h 0.03017 cos 46» - 0.00476 cos 661 -h 



(A23) 



where it is sufficient to consider the first six terms of Eq. 
(A23) for an accuracy of < 10~^. It is important to note 
that Eq. (A21) merely constitutes a special solution to 
Eq. (A5), and thus we still need to address the question 
of boundary conditions in the deep-MOND case. 

Given a suitable choice of such boundary conditions, 
Eq. (A21) needs to be replaced by 

$QM (r, = Ar^l'^T{e) + $0 (r, + const , ( A24) 

where ^^{r.Q) is an appropriate solution of Laplace's 
equation. Using the previously assumed smoothness and 
periodicity condition of the angular part and additionally 
requiring that the potential is regular at r = 0, the most 
general form of '^air, 9) can be expressed as 



n=l 



2n-l 



r2"P2„(c0S^), 



n e 



(A25) 



where C is defined according to Eq. (A2), P„ denotes the 
Legendre polynomial of degree n and g„ € M. Since the 
validity of equation (A24) is restricted to a close vicinity 
of the SP at r = 0, however, wc expect that only the first 
few terms in Eq. (A25) will significantly contribute to 
the solution in this case. As we shall see later, this allows 
one to approximate the infinite series in Eq. (A25) by a 
cut-off one. Unfortunately, there exists no independent 
way of determining the coefficients (?„ because the saddle 
region is (i(n'oid of any real sources. Thus their values 
are set by the boundary conditions in the intermediate 
MOND domain and need to be estimated from the full 
numerical result, which is further discussed in the next 
section. 



3. Comparison to numerical results 

Having derived semi-analytic solutions in both the 
quasi-Newtonian and the deep-MOND domain, we now 

want to compare these to the full numerical two-body 
results of the Earth-Sun SP. To do so, however, we still 



le-15 - 



X [km] 



FIG. 18. Potential derivative Syy of the Earth-Sun SP eval- 
uated at J/ « 1.4 and z — in units of km: Shown are the 
numerical result for u{y) given by Eq. (12) with k = 0.03 
(solid line) as well as the semi-analytic approximations in the 
quasi-Newtonian (dashed line) and deep-MOND (dotted line) 
regime, respectively. 



need to determine the series coefficients in Eq. (A6) for a 
particular choice of J^(y). Let us consider the expression 
in Eq. (12) which is the QMOND analog of the TeVeS 
interpolating fimction used in Ref. [18]. For y ^ 1, the 
above formula can be expanded as 



1 /47r 

4 It 



yZ 



O 



y" 



(A26) 



which reveals that the first-order term vanishes, i.e. one 

has ai = 0. Since the numerical approach uses different 
boundary conditions (see Sec. IV B), we choose not to 
compare the actual potentials, but rather their second 
derivatives. In accordance with our numerical setup, wc 
switch to Cartesian coordinates, letting the symmetry 
axis now point along the a;-direction. 

As for a comparison to the deep-MOND solution, wc 
must also take care of the homogeneous contribution $0 
which is determined by the potential values in the in- 
termediate MOND domain. The coefficients g„ can be 
calculated by fitting the simulation results with the help 
of Eq. (A24), which, of course, requires introducing a 
cut-off in Eq. (A25). Since only the first few terms 
signific:antly contribute in the close vicinity of the SP, 
wc cut the scries off at n = 3 and discard all remain- 
ing terms. Furthermore, to ensure that Eq. (A18) is 
approximately satisfied, we consider only data points 
with 0.05 < r/ro < 0.5 for this procedure [18]. Set- 
ting k = 0.03 in Eq. (12) and using the corresponding 
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FIG. 19. Same as Fig. 18, but now for the function vly) given 
by Eq. (11). 



simulation data of the Earth-Sun SP, we find 



qi RJ 2.001 X 10-^ 

92 « -2.433 X 10-^^ 

93 w 1.228 X 10"^^. 



(A27) 
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FIG. 20. Rotation curves of NGC 1560 (upper panel), 
NGG2403 (middle panel) and NGG 4157 (lower panel): hile 
The data points correspond to observations. The dashed, 
solid and dotted lines represent rough "fits" obtained using 
Eq. (12) , Eq. (11) and Eq. (13) and setting MjL = 1 
(for NGC 1560 and NGC 2403) and MjL = 1.5 (for NGC 
4157), respectively. Also shown are the Newtonian rotation 
curves of the individual stellar (dot-dashed line) and gas (long 
dot-dashed line) components. 



Fixing y « 1.4 and z = in units of km, the numeri- 
cal two-body result for Syy = dy<^QM of the Earth-Sun 
SP is illustrated in Fig. 18. Since the numerical result 
is essentially symmetric with respect to x = 0, we only 
consider the case a; > 0. As can be seen from the figures, 
the approximate analytic solutions are well recovered by 
the full numerical result. Note that the increasing devia- 
tion from the quasi-Newtonian solution near the bound- 
ary at a; = lO^km is entirely due to the different choice 
of boundary conditions (see Sec. IV). The numerical 
solution becomes unreliable for x < 20km, which is a 
consequence of the grid's finite resolution. 

Another example is given by the simple interpolating 
function defined in Eq. (11). Its expansion takes the 
form 

Prom the above, one immediately identifies ai = 1 and 
a2 = — 1- As is shown in Fig. 19 for Syy, the quasi- 
Newtonian solution provides an excellent approximation 
to the full numerical result in this case. On the other 
hand, the deep-MOND limit is not reached at the given 
resolution (Axmax ~ 2.8km). 



Appendix B: Impact of 1^(3/) on galactic scales 

Any functional form of i^{y) chosen for the SS is re- 
quired to be consistent with the vast amount of available 
extragalactic data. The properties of Eqs. (11) and (12) 
relevant to these scales are well-known and have been ex- 
tensively discussed in the literature (see, e.g., Ref. [45]). 
While both prescriptions are able to describe the ob- 
served rotation curves of spiral galaxies, it is usually con- 
cluded that the data prefers the simple choice Eq. (11) 
(or further refined versions thereof) since it provides a 
better fit to observations if one leaves the stellar mass- 
to-light ratio M/L as a free parameter. For Eq. (13), 
however, the situation is unclear. To get an idea about its 
performance, we consider the rotation curves of three test 
galaxies covering different gravitational regimes: NGC 
1560 (low surface brightness), NGC 2403 and NGC 4157 
(high surface brightness). For the different choices of 
v{y), Fig. 20 illustrates the predicted MOND rotation 
curves together with the observed ones, assuming con- 
stant M/L ratios in the B-band. Note that since we are 
only interested in the performance of Eq. (13) compared 
to the others, we have not tried to actually fit the ob- 
served curves. Instead, to obtain a very rough match to 
the data, we have used M/L = 1 for both NGC 1560 and 
NGC 2403 and M/L = 1.5 for NGC 4157. As can be seen 
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FIG. 21. Formally constructed TeVeS interpolating functions 
ji for k = 0.03 corresponding to the i/-functions given by Eq. 
(12) (solid line), Eq. (11) (dashed line) and Eq. (13) (dotted 
line) . 



from the figure, the behavior of Eq. (13) is very similar 
to that of the simple j/-function, despite their differences 
on SS scales. This suggests that the interpolating func- 
tion Eq. (13) is not only compatible with extragalactic 
data, but also recovers the successful descriptive power 
ofEq. (11). 



Appendix C: Interpolating functions for TeVeS 

So far we have only considered the situation of low- 
acceleration bubbles within the QMOND framework 
which emerges in the nonrelativistic limit of particular 
bimetric theories (see Sec. II). Another possibility of 
recovering MONDian dynamics in this limit is realized 
in TeVeS. In the following, we want to address whether 
there exist viable interpolation functions for the TeVeS 
scalar field 4> which give rise to similar results as found 
in Sec. IV. 

Applying the approximations for weak fields and quasi- 
static systems, the gradient of the total gravitational po- 
tential in TeVeS approximately takes the form [14, 46] 



V<I> = V<I> 



N 



(CI) 



where $Ar is the Newtonian potential, obeying the usual 
Poisson equation, and the scalar field (/> is a solution of 



kGp. 



(C2) 



Here p denotes the density of real matter and fi depends 
on the scalar field strength gs = |V(/)| [51]. In situations 



where the curl field vanishes, it follows from Eq. (C2) 
that 



k 



-V4> 



N 



which, if combined with Eq. (CI), leads to 



V$ = 1 



N- 



(C3) 



(C4) 



Comparing the above with Eq. (1), we identify the 
QMOND interpolating function v as 



k 



(C5) 



Together with our previous assumption that v ap- 
proaches unity for strong gravitational fields, the above 
leads to /x — > oo in this limit. Such interpolation func- 
tions jl have been found to be problematic and may there- 
fore be regarded as inapplicable [52] . Avoiding this issue 
by imposing /i — > 1 for g^/ao ^ 1, we obviously must 
have that v ^ 1 + k/ [Air) in the Newtonian limit, consis- 
tent with Eq. (C5) and resulting in an effective rescaling 
of the gravitational constant. Since we are here interested 
in the purely non-Newtonian part of the gravitational 
potential, this can be easily achieved by formally adding 
the constant A:/(47r) to a given expression for v{y). As 
a consequence, the full solutions of Eq. (1) for the such 
constructed i/- functions only differ by the term k^N / (47r) 
from their counterparts, which leads to the same results 
as in Sec. IV after subtracting the total Newtonian con- 
tribution. Another way of saying this is the following: 
Adding a constant to v{y) does not alter the effective 
density field if the region of interest is devoid of any real 
physical matter sources. Choosing the same boundary 
conditions for the potential then yields identical results. 
Keeping this technical detail in mind, we are now ready 
to construct the corresponding functions jl with the help 
of Eqs. (C3) and (C5). Starting with the interpolation 
function Eq. (12), a bit of algebra reveals that 



k 



yj\ - jl'^ 47rao 



9s 



(C6) 



which is exactly the implicit definition of /i used in Refs. 
[16, 18]. Similarly, one may find expressions for the in- 
terpolating functions given by Eqs. (11) and (13). As- 
suming k ~ 0.03, Fig. 21 shows the resulting TeVeS 
interpolation functions /i for 0.1 < gs/da S lO"'- While 
all of them exhibit the same asymptotic behavior, the re- 
sulting /i((7s/ao) for Eq. (13) (dotted line) is multivalued, 
and therefore ill-defined [45, 52]. Since it is unclear how 
to assign /t-values to a given g^ in this case, the boundary 
value problem in Eq. (C2) turns ambiguous, rendering 
the theory unphysical. The only remedy is to resort to a 
larger renormalization of the gravitational constant, i.e. 
a larger fc, but this would result in unacceptable cos- 
mological predictions [28]. The interpolating functions 
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constructed from Eqs. (11) and (12) (dashed and solid 
lines) do not suffer from this problem and lead to healthy 
theories. The corresponding TeVeS stress signals will be 
comparable to the QMOND results, but slightly lower 



due to a softening caused by the curl field [28], which 
allows one to adopt the conclusions of Sec. IV for these 
cases. 
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